3.Image Analysis

library(data.table)
library(anytime)
library(stringr)
library(dplyr)
library(gridExtra)
library(ggplot2)
library(plotly)
library(jpeg)
library(fields)
library(formatR)

1. Read Image

img = readJPEG("/Users/ilkerkurtulus/Downloads/IMG_7955.jpg")

str(img)

##  num [1:512, 1:512, 1:3] 0.424 0.4 0.388 0.404 0.4 ...

dim(img)

## [1] 512 512   3

if (exists("rasterImage")) {
    plot(1:2, type = "n")
    rasterImage(img, 1, 1, 2, 2)
}

if (exists("rasterImage")) {
    plot(1:4, type = "n")
    rasterImage(img[, , 1], 1, 1, 2, 2)
    rasterImage(img[, , 2], 2, 1, 3, 2)
    rasterImage(img[, , 3], 3, 1, 4, 2)
}

noise = runif(512 * 512 * 3, min = 0, max = 1)
dim(noise) = c(512, 512, 3)
nimg = noise + img
nimg_scaled = (nimg - min(nimg))/(max(nimg) - min(nimg))

if (exists("rasterImage")) {
    plot(1:2, type = "n")
    rasterImage(nimg_scaled, 1, 1, 2, 2)
}

if (exists("rasterImage")) {
    plot(1:4, type = "n")
    rasterImage(nimg_scaled[, , 1], 1, 1, 2, 2)
    rasterImage(nimg_scaled[, , 2], 2, 1, 3, 2)
    rasterImage(nimg_scaled[, , 3], 3, 1, 4, 2)
}

4. Define greyscale function according to luminosity type:

to_greyscale = function(img) {
    img[, , 1] * 0.21 + img[, , 2] * 0.72 + img[, , 3] * 0.07
}
grey_nimg = to_greyscale(nimg_scaled)

pca_img = prcomp(grey_nimg, center = TRUE, scale. = TRUE)
eigs = pca_img$sdev^2
exp_var_ratio = eigs/sum(eigs)
cum_exp_var_ratio = cumsum(exp_var_ratio)
plot(cum_exp_var_ratio)

Even the image has 512 dimensions we can explain it well with 250 components with 0.95 confidence interval. However its still very large dimension for statistical inference. b) Construct patch matrix:

extract_path = function(data, i, j, n) {
    dw = (n - 1)/2
    as.vector(t(data[(i - dw):(i + dw), (j - dw):(j + dw)]))
}
tmp = c()
for (i in 2:511) {
    for (j in 2:511) {
        x = extract_path(grey_nimg, i, j, 3)
        tmp = c(tmp, x)
    }
}
res = matrix(tmp, nrow = 260100, byrow = TRUE)

Apply pca and plot:

pca_res = prcomp(res, center = TRUE, scale. = TRUE)

# 1st component
pca_res1d = pca_res$x[, 1]
pca_mat = matrix(pca_res1d, nrow = 510, ncol = 510, byrow = TRUE)
pca_mat_scaled = (pca_mat - min(pca_mat))/(max(pca_mat) - min(pca_mat))
if (exists("rasterImage")) {
    plot(1:2, type = "n")
    rasterImage(pca_mat_scaled, 1, 1, 2, 2)
}

# 2nd component
pca_res1d = pca_res$x[, 2]
pca_mat = matrix(pca_res1d, nrow = 510, ncol = 510, byrow = TRUE)
pca_mat_scaled = (pca_mat - min(pca_mat))/(max(pca_mat) - min(pca_mat))
if (exists("rasterImage")) {
    plot(1:2, type = "n")
    rasterImage(pca_mat_scaled, 1, 1, 2, 2)
}

# 3rd component
pca_res1d = pca_res$x[, 3]
pca_mat = matrix(pca_res1d, nrow = 510, ncol = 510, byrow = TRUE)
pca_mat_scaled = (pca_mat - min(pca_mat))/(max(pca_mat) - min(pca_mat))
if (exists("rasterImage")) {
    plot(1:2, type = "n")
    rasterImage(pca_mat_scaled, 1, 1, 2, 2)
}

ev1 = pca_res$rotation[, 1]
ev_scaled1 = (ev1 - min(ev1))/(max(ev1) - min(ev1))
mat_ev1 = matrix(ev_scaled1, nrow = 3, ncol = 3, byrow = TRUE)

ev2 = pca_res$rotation[, 2]
ev_scaled2 = (ev2 - min(ev2))/(max(ev2) - min(ev2))
mat_ev2 = matrix(ev_scaled2, nrow = 3, ncol = 3, byrow = TRUE)

ev3 = pca_res$rotation[, 3]
ev_scaled3 = (ev3 - min(ev3))/(max(ev3) - min(ev3))
mat_ev3 = matrix(ev_scaled3, nrow = 3, ncol = 3, byrow = TRUE)


if (exists("rasterImage")) {
    plot(1:4, type = "n")
    rasterImage(mat_ev1, 1, 1, 2, 2)
    rasterImage(mat_ev2, 2, 1, 3, 2)
    rasterImage(mat_ev3, 3, 1, 4, 2)
}